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RUNOFF  SIMULATION  USING  RADAR  RAINFALL  DATAi 


John  C.  Peters  and  Daniel  J.  Easton‘S 


ABSTRACT:  Rainfall  data  products  generated  with  the  national 
network  of  WSR-88D  radars  are  an  important  new  data  source  pro¬ 
vided  by  the  National  Weather  Service.  Radar-based  data  include 
rainfall  depth  on  an  hourly  basis  for  grid  cells  that  are  nominally  4 
km  square.  The  availability  of  such  data  enables  application  of 
improved  techniques  for  rainfall-runoff  simulation.  A  simple  quasi - 
distributed  approach  that  applies  a  linear  runoff  transform  to  grid- 
ded  rainfall  excess  has  been  developed.  The  approach  is  an 
adaptation  of  the  Clark  conceptual  runoff  model,  which  employs 
translation  and  linear  storage.  Data  development  for,  and  results 
of,  an  initial  application  to  a  4160  km^  watershed  in  the  Midwest¬ 
ern  U.S.  are  illustrated. 

(KEY  TERMS:  hydrograph  analysis  and  modeling;  simulation;  sur¬ 
face  water  hydrology;  radar.) 


INTRODUCTION 

Traditional  application  of  the  unit  hydrograph 
approach  to  runoff  simulation  involves  the  use  of  spa¬ 
tially  averaged  (lumped)  values  of  basin  rainfall  and 
infiltration  (losses).  This  approach  has  been  of  practi¬ 
cal  value  because  data  available  from  typically  sparse 
rain-gage  networks  are  generally  inadequate  to  justi¬ 
fy  more  spatially  detailed  simulation  methods.  The 
availability  of  “new-generation”  radar  rainfall  data 
enhances  the  attractiveness  of  distributed  simulation 
approaches  that  take  into  account  spatial  variations 
of  rainfall  and  watershed  characteristics. 

To  facilitate  initial  use  of  radar  rainfall  data,  a  rel¬ 
atively  simple  quasi-distributed  approach  has  been 
developed  that  applies  a  linear  runoff  transform  to 
gridded  rainfall  excess.  The  approach  is  an  adaptation 
of  the  Clark  conceptual  runoff  model  (Clark.  1943), 
which  represents  surface  runoff  with  translation  and 
linear-storage  attenuation.  In  this  adaptation,  radar 


grid  cells  are  superposed  on  the  basin,  and  rainfall 
and  losses  are  tracked  uniquely  for  each  cell.  Rainfall 
excess  for  each  cell  is  lagged  to  the  basin  outlet  by  the 
cell’s  travel  time  (i.e.,  time  of  travel  from  the  cell  to 
the  basin  outlet).  The  lagged  excesses  are  routed 
through  a  linear  reservoir,  and  baseflow  is  added  to 
obtain  a  total-runoff  hydrograph.  The  computer  pro¬ 
gram  that  performs  these  computations  is  the  Modi¬ 
fied  Clark  (modClark)  Runoff  Simulation  Program 
(HEC,  1995a). 


RADAR  RAINFALL  DATA 

A  national  network  of  WSR-88D  (Weather  Surveil¬ 
lance  Radar-88  Doppler)  radars  is  being  deployed  by 
the  National  Weather  Service  (NWS).  Processing  of 
precipitation  data  by  the  NWS  is  done  in  stages 
(Shedd  and  Fulton,  1993).  Stage  III  products  incorpo¬ 
rate  information  from  “ground  truth”  rain  gages  and 
satellite  and  surface  temperature  observations,  and 
they  result  from  merging  (“mosaicking”)  data  from 
overlapping  radar  coverages.  For  the  application  illus¬ 
trated  subsequently.  Stage  III  hourly  precipitation 
data  were  obtained  via  Internet  from  the  NWS 
Arkansas-Red  Basin  River  Forecast  Center  (ABRFC) 
in  Tulsa,  Oklahoma.  As  of  mid-1995,  the  ABRFC  was 
the  only  NWS  River  Forecast  Center  from  which 
Stage  III  products  were  routinely  available,  although 
testing  of  Stage  III  processing  was  underway  at  sever¬ 
al  other  River  Forecast  Centers. 

The  Stage  III  rainfall  data  are  provided  for  cells 
defined  by  the  Hydrologic  Rainfall  Analysis  Project 
(HRAP)  grid  (Greene  and  Hudlow,  1982).  The  HRAP 


^Paper  No.  95150  of  the  Water  Resources  Bulletin.  Discussions  are  open  until  February  1, 1997. 

Respectively.  Senior  Engineer,  U.S.  Army  Corps  of  Engineers,  Hydrologic  Engineering  Center,  609  Second  St.,  Davis,  California  95616- 
4687;  and  Graduate  Student,  Department  of  Civil  Engineering,  University  of  California,  Davis,  California  95616. 


1 


Peters  and  Easton 


gnd  is  uniform  on  a  polar  stereographic  map  projec- 
tion.  Consequently,  the  dimensions  of  an  HRAP  grid 
cell,  as  projected  on  the  earth’s  surface,  vary  with  lati¬ 
tude.  Figure  1  illustrates  an  HRAP  grid  superposed 
over  four  subbasins  of  the  4160  km2  Illinois  River 
watershed  upstream  from  Tenkiller  Lake.  The  water¬ 
shed  is  located  in  northeastern  Oklahoma  and  north¬ 
western  Arkansas.  The  grid  cell  areas  vary  in  this 
watershed  from  16.3  to  16.5  km^. 


Figure  1.  HRAP  Grid  Superposed  on  Four 
Subbasins  of  the  Illinois  River  Watershed. 


Radar  rainfall  data  obtained  from  the  ABRFC  is  in 
the  netCDF  (Network  Common  Data  Form)  format 
(Unidata  Progp’am  Center,  1991).  A  utility  program 
titled  gridUtl  (HEC,  1995b)  loads  the  data  into  a 
direct  access  file  associated  with  the  Hydrologic  Engi¬ 
neering  Center’s  Data  Storage  System  (HEC-DSS). 
The  Modified  Clark  program  retrieves  the  gridded 
rainfall  data  from  an  HEC-DSS  file. 


MODIFIED  CLARK  METHOD 

'Two  basin  parameters  are  required  to  transform 
rainfall  excess  to  direct  runoff  with  the  Modified 
Clark  method:  time  of  concentration,  T^-,  and  storage 
coefficient  (for  a  linear  reservoir),  R.  Both  have  units 
of  time.  'Translation  is  performed  on  a  grid  cell  basis 
by  using  a  travel  time  index.  The  travel  time  (or 
translation  lag)  for  a  grid  cell  is  calculated  as  follows: 

(.travel  =  T,  (1) 

{travel  time  index)nax 

where  is  the  time  of  concentration  for  the  basin, 
{travel  time  index^^n  is  the  travel  time  index  for  a 
cell,  and  {travel  time  index),^^  is  the  maximum  travel 
time  index  of  all  of  the  cells  associated  with  the  basin. 
The  development  of  a  travel  time  index  is  described  in 
the  next  section. 

The  lagged  rainfall  excess  for  each  cell  is  routed 
through  a  linear  reservoir  with  the  following  equa¬ 
tion: 


T  -L 

1  _ 

_i2  +  0.5*A^  J 

^avg  ^ 

JL  — 

i?  +  0. 5* 

where  Oi  is  direct  runoff  at  time  i,  R  is  the  storage 
coefficient,  is  the  average  inflow  for  the  interval 
i-1  to  i,  and  Af  is  the  time  interval. 


CELL  PARAMETERS 

Part  of  the  required  input  for  the  Modified  Clark 
program  is  a  cell-parameter  file  that  contains  the 
following  information  for  each  cell:  cell  x-coordinate, 
cell  y-coordinate,  area  (within  basin),  and  travel  time 
index.  As  shown  in  the  previous  section,  the  travel 
time  index  for  a  cell  is  used  to  calculate  a  translation 
lag.  The  travel  time  from  a  cell  to  the  basin  outlet  is 


where  x  is  the  time-of-travel  to  the  basin  outlet,  D  is 
the  length  of  the  flow  path  to  the  basin  outlet,  and 
Va„g  is  the  average  velocity  over  the  flow  path.  If  it  is 
assumed  that  travel  velocity  is  constant  for  the  basin, 
then  flow  path  length  can  serve  as  the  cell  travel  time 
index. 

An  alternative  to  the  assumption  of  a  constant 
travel  velocity  is  to  incorporate  a  spatially  distributed 
velocity  field,  as  proposed  by  Maidment  et  al.  (1996). 
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The  travel  velocity  through  a  cell  is  assumed  to  be 
proportional  to  the  cell  slope  and  to  the  accumulated 
area  of  all  cells  contributing  runoff  to  the  cell.  That  is, 

^ceil  (4) 

where  Vceii  is  the  travel  velocity  through  a  cell,  S  is 
the  cell  slope,  and  A  is  the  accumulated  area  of  con¬ 
tributing  cells.  The  accumulated  area  can  be  regarded 
as  a  surrogate  for  depth.  A  value  of  0.5  has  been 
found  to  be  reasonable  for  both  the  a  and  b  exponents 
(Maidment  et  ai,  1996).  The  travel  time  index  for  a 
cell  is  then  defined  as  the  integral  of  ^ceZ/^^cell  along 
the  flow  path  to  the  basin  outlet,  where  Iceii  is  the 
length  of  flow  path  through  a  cell.  Incorporation  of  a 
spatially  distributed  velocity  field  in  computing  travel 
time  indices  is  worthy  of  further  study.  However,  for 
the  purposes  of  this  paper,  the  assumption  of  a  con¬ 
stant  average  velocity  over  all  the  basin  flow  paths  is 
adopted  for  an  initial  demonstration  of  the  Modified 
Clark  method. 

Procedures  for  using  a  geographic  information  sys¬ 
tem  (GIS)  to  calculate  cell  areas  and  travel  time 
indices  have  been  developed  (HEC,  1995c).  The  proce¬ 
dures  require  processing  digital  elevation  model 
(DEM)  data  such  as  are  available  for  the  continental 
U.S.  (via  Internet)  from  the  USGS  EROS  Data  Center 
(USGS,  1990).  An  eight-direction  “pour-point”  algo¬ 
rithm  defines  the  direction  of  flow  from  any  grid  cell 
to  be  in  the  direction  of  steepest  descent  from  the  cell 
to  one  of  its  eight  neighbors.  A  flow  path  length  is 
computed  by  summing  the  lengths  of  all  segments 
along  the  path  from  the  cell  to  the  basin  outlet.  Area 
and  travel  time  index  are  determined  for  DEM-based 
cells  at  a  100  m  resolution.  Radar  cells  (based  on  the 
HEAP  grid)  are  then  superposed  and  their  areas  and 
travel  time  indices  are  calculated  by  summing  the 
areas  and  averaging  the  travel  time  indices  of  the 
encompassed  DEM-based  cells.  The  cell  areas  and 
travel  time  indices  are  treated  as  constants  for  a 
given  basin.  Thus  GIS  is  used  for  a  one-time  process¬ 
ing  of  data  and  is  not  required  for  subsequent  applica¬ 
tion  of  the  Modified  Clark  program. 


LOSSES,  BASEFLOW,  AND 
HYDROLOGIC  ROUTING 


Loss  models  available  in  the  Modified  Clark  pro¬ 
gram  are  Initial/Constant,  SCS  Curve  Number,  and 
Green  and  Ampt.  The  methods  are  applied  as  in  the 
HEC-1  program  (HEC,  1990).  The  loss  model  parame¬ 
ters  apply  to  all  cells  in  the  basin,  but  losses  are 
calculated  individually  for  each  cell  based  on  the  rain¬ 
fall  intensities  associated  with  that  cell.  Baseflow  is 


modeled  as  in  HEC-1,  The  starting  flow,  recession 
flow,  and  recession  ratio  parameters  are  used  to  calcu¬ 
late  baseflow  at  the  outlet  of  the  basin. 

The  Modified  Clark  program  can  only  simulate 
runoff  from  elemental  basins  -  that  is,  basins  that  are 
not  subdivided.  However,  the  program  has  the  capa¬ 
bility  to  write  its  simulation  results  (i.e.,  discharge 
hydrographs)  to  the  HEC-DSS.  For  applications  with 
multi-subbasin  watersheds,  the  hydrographs  can  be 
retrieved  from  HEC-DSS  and  routing  performed  with 
programs  such  as  HEC-1,  HEClF  (Peters  and  Ely, 
1985),  or  UNET  (HEC,  1993). 


TEST  WATERSHED 

Runoff  simulations  were  performed  for  the  Illinois 
River  watershed  above  Tenkiller  Lake  in  northeastern 
Oklahoma  and  northwestern  Arkansas.  The  4,163 
km2  watershed  was  divided  into  four  subbasins 
as  shown  in  Figure  2.  The  subbasin  areas  and  the 
number  of  radar  cells  in  each  subbasin  are  listed  in 
Table  1.  Stream  gages  are  located  at  the  outlets  of 
subbasins  85,  86,  and  113.  Inflow  to  Tenkiller  Lake 
can  be  computed  from  measured  outflow  and  lake- 
level  data.  Figure  2  also  shows  the  location  of  precipi¬ 
tation  gages,  for  which  hourly  rainfall  is  available. 
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Figure  2.  Illinois  River  Watershed. 
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TABLE  1.  Subbasin  Area  and  Number  of  Radar  Cells. 


Subbasin 

Area 

Ckm2) 

Number  of 
Radar  CeUs 

85 

829 

84 

86 

1645 

129 

113 

795 

78 

127 

894 

86 

The  watershed  is  in  the  Ozark  Highlands  and  is 
heavily  wooded.  Elevations  range  from  140  meters 
above  sea  level  at  the  outlet  of  Tbnkiller  Lake  to  580 
meters.  The  hills  in  the  region  are  formed  of  porous 
limestone  and  overlain  with  cherty  topsoil.  The  flood 
plains  can  be  gravelly,  and  in  places  the  substratum  is 
too  pervious  to  hold  water.  Therefore,  high  infiltration 
is  expected  (Soil  Conservation  Service,  1965  and 
1970).  For  simplicity,  the  method  of  using  an  initial 
loss  followed  by  a  constant  loss  rate  was  adopted  for 
calculating  rainfall  excess. 


STORM  EVENTS 

Radar  rainfall  data  for  storms  that  occurred  on 
November  4-5,  1994,  January  13-14, 1995,  and  May  8, 
1995,  were  used  for  the  initial  application  of  the  Mod¬ 
ified  Clark  method.  Table  2  shows  total  average  rain¬ 
fall  for  each  storm  over  the  four  subbasins  as 
calculated  using  (a)  Stage  III  radar  data  and  (b)  data 
from  the  precipitation  gages  shown  in  Figure  2.  Total 
average  rainfall  from  the  gage  data  was  calculated 
^sing  an  inverse  distance-SQuared  weighting  proce¬ 
dure  (HEC,  1989). 

The  total  average  precipitation  calculated  for  each 
of  the  three  storm  events  using  gage  data  differs  sig¬ 
nificantly  from  that  calculated  using  radar  data.  Dif¬ 
ferences  might  be  attributed  to  various  factors, 
including  the  spatial  variability  of  the  rainfall* 
weighting  of  the  gage  data,  the  accuracy  of  the  radar 
rainfall  data,  and  associated  processing  procedures. 
While  these  are  key  issues  with  regard  to  rainfall 
measurement,  their  resolution  is  beyond  the  scope  of 
this  paper,  which  is  intended  to  demonstrate  use  of 
the  gridded  rainfall  data. 

A  time-area  concentration  histogram  for  subbasin 
85  is  shown  in  Figure  3.  The  histogram  is  based  on 
the  area  and  travel  time  index  for  each  radar  cell,  and 
it  shows  the  percent  of  the  subbasin  area  that 
contributes  runoff  at  the  outlet  (via  translation)  for 
increments  of  travel  time  (expressed  as  10  percent 
increments  of  the  time  of  concentration).  A  time- 
volume  concentration  histogram  for  the  November 


4-5,  1994,  storm,  which  shows  the  percent  of  the  total 
volume  of  rainfall  that  contributes  runoff  to  the  outlet 
for  increments  of  travel  time,  is  also  shown  in  Figure 
3.  If  the  rainfall  were  distributed  uniformly,  the  two 
histograms  would  be  identical;  the  histograms  differ 
because  of  spatial  variations  in  rainfall.  As  shown, 
the  time-volume  histogram  does  not  vary  greatly  from 
the  time-area  histogram.  This  was  generally  true  for 
the  three  storm  events  over  the  Illinois  River  water¬ 
shed.  Conclusions  about  the  time-space  rainfall  distri¬ 
bution  cannot  be  made  from  these  histograms  because 
the  hourly  cell  data  have  been  integrated  over  time. 


TABLE  2.  Total  Average  Rainfall  as  Calculated 
Using  Radar  and  Gage  Rainfall  Data. 


Storm  Event 

Subbasin 

Total 

Average 

Rainfall 

(radar) 

(mm) 

Total 

Average 

Rainfall 

(gage) 

(mm) 

November  4-5,  1994 

85 

93 

74 

November.  4-5,  1994 

86 

90 

November  4-5,  1994 

113 

98 

118 

November  4-5,  1994 

127 

92 

133 

Januaiy  13-14,  1995 

85 

91 

51 

January  13-14,  1995 

86 

94 

44 

January  13-14,  1995 

113 

81 

55 

January  13-14, 1995 

127 

66 

75 

May  8,  1995 

85 

57 

37 

May  8,  1995 

86 

56 

31 

May  8,  1995 

113 

58 

39 

May  8, 1995 

127 

62 

81 

%  of  Time  of  Concentration 


Figure  3.  Time-.Area  and  Time-Rainfall  Volume  Histograms 
for  the  November  4-5,  1994,  Storm  on  Subbasin  85. 
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MODIFIED  CLARK  SIMULATION 

Results  of  the  Modified  Clark  runoff  simulations  at 
the  Watts,  Tahlequah,  and  Eldon  gages  and  Tenkiller 
Lake  for  the  November  4-5,  1994,  January  13-14, 
1995,  and  May  8,  1995,  storms  are  shown  in  Figures 
4,  5,  and  6,  respectively.  Loss  parameters  were  adjust¬ 
ed  so  that  the  volumes  of  observed  and  simulated 
runoff  were  essentially  identical.  The  Clark  (i.e., 
basin  time-of-concentration  and  storage  coefficient), 
loss,  and  baseflow  parameters  used  in  the  simulations 
are  shown  in  Table  3.  Values  for  time-of-concentration 
and  storage  coefficient  were  kept  constant  for  the  sim¬ 
ulations.  Flow  simulation  at  the  Tahlequah  gage  sta¬ 
tion  and  Tenkiller  Lake  required  stream  routing  of 
hydrographs  generated  at  upstream  locations.  This 
was  performed  using  the  modified  Puls  method  as 
implemented  in  HEC-1  (HEC,  1990)  with  storage-dis¬ 
charge  criteria  furnished  by  the  Tulsa  District  of  the 
Corps  of  Engineers. 


As  shown  in  Figures  4,  5,  and  6,  the  simulated 
hydrographs  provide  a  reasonable  fit  to  the  observed 
hydrographs.  Simulations  were  also  performed  using 
spatially  averaged  radar-rainfall  data.  The  results 
were  similar  to  those  based  on  grid-distributed  rain¬ 
fall.  This  is  attributed  to  the  uniformity  of  the  rainfall 
distribution  as  discussed  in  the  previous  section.  It  is 
expected  that  with  an  application  to  a  storm  with 
marked  spatial  variability,  such  as  a  localized  convec¬ 
tive  storm,  a  substantial  difference  would  occur 
between  simulations  based  on  grid-distributed  versus 
spatially-averaged  rainfall.  The  difference  would  be 
due  to  both  the  grid-based  calculation  of  losses  as  well 
as  the  grid-based  translation  of  rainfall  excess.  Hypo¬ 
thetical  data  have  been  used  to  confirm  this  conclu-  ^ 

sion,  but  data  have  not  been  available  for  the  ^ 

watershed  above  Tenkiller  Lake  for  such  comparisons. 
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Figure  4.  Modified  Clark  Rainfall-Runofi*  Simulations  for  the  November  4-5,  1994,  Storm. 
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Figure  5.  Modified  Clark  Rainfall-Runoff  Simulations  for  the  January  13-14,  1995,  Storm. 


CONCLUDING  REMARKS 

The  availability  of  rainfall  data  from  WSR-88D 
radars  affords  new  opportunities  for  increasing  the 
spatial  detail  with  which  rainfall-runoff  processes  are 
.  simulated.  A  simple  method  for  simulating  watershed 
runoff  by  using  a  linear  transform  of  grid-distributed 
rainfall  excess  is  described  herein.  Aside  from  cell 
properties  (which  can  be  obtained  with  GIS  proce¬ 
dures),  the  data  requirements  for  the  Modified  Clark 
method  are  essentially  the  same  as  for  existing 
lumped-parameter  models.  The  method  thus  provides 
a  relatively  straightforward  transition  to  use  of  radar- 
rainfall  data.  As  more  physically  based  distributed 
models  come  into  use,  it  may  be  useful  to  compare 
their  performance,  data  requirements,  and  utility 
with  a  simpler  approach  such  as  that  described 
herein. 
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Figure  6.  Modified  Clark  Rainfall-Runoff  Simulations  for  the  May  8,  1995,  Storm. 


TABLE  3.  Modified  Clark  Model  Parameters  Used  in  Test  Simulations, 


Subbasin 

Storm 

Event 

Clark  Parameters 

Loss  Parameters 

Baseflow  Parameters 

Tc 

Giours) 

R 

(hours) 

Initial 

Loss 

(mm) 

Constant 
Loss  Rate 
(nun/hr) 

Initial 
Flow 
(cu.  m/s) 

Recession 

Ratio* 

85 

November  4-5,  1994 

'  30 

15.5 

62.3 

1.3 

1.4 

1.02 

January  13-14,  1995 

30 

15.5 

57.2 

1.0 

2.8 

1.02 

May  8,  1995 

30 

15.5 

17.5 

0,0 

0.0 

1.02 

86 

November  4-5,  1994 

24 

11.6 

69.9 

3.6 

8.5 

1.02 

January  13-14,  1995 

24 

11.6 

27.9 

5.1 

8.5 

1.02 

May  8,  1995 

24 

11.6 

30.2 

0.0 

113.0 

1.02 

113 

November  4-5,  1994 

16 

9.7 

92.7 

1.8 

1.4 

1.02 

January  13-14,  1995 

16 

9.7 

56.6 

0.0 

2.8 

1.02 

May  8,  1995 

16 

9.7 

21.1 

0.0 

65.1 

1.02 

127 

November  4-5,  1994 

1 

10.0 

32.0 

1.5 

0.0 

1,00 

January  13-14,  1995 

1 

10.0 

30.0 

0.0 

0.0 

1.00 

May  8,  1995 

1 

10.0 

12.7 

0.0 

0.0 

1.00 

*The  ratio  is  that  of  the  initial  flow  to  the  flow  one  hour  later. 
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